knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "~/DADS")

library(phyloseq); packageVersion("phyloseq")
## [1] '1.50.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.5.1'
library(readr); packageVersion("readr")
## [1] '2.1.5'
library(purrr); packageVersion("purrr")
## [1] '1.0.4'
library(furrr); packageVersion("furrr")
## Loading required package: future
## [1] '0.3.1'
library(dplyr); packageVersion("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## [1] '1.1.4'
library(stringr); packageVersion("stringr")
## [1] '1.5.1'
library(metacoder); packageVersion("metacoder")
## This is metacoder version 0.3.8 (stable)
## 
## Attaching package: 'metacoder'
## The following object is masked from 'package:ggplot2':
## 
##     map_data
## The following object is masked from 'package:phyloseq':
## 
##     filter_taxa
## [1] '0.3.8'
library(data.table); packageVersion("data.table")
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:purrr':
## 
##     transpose
## [1] '1.16.4'
library(Biostrings); packageVersion("Biostrings")
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:data.table':
## 
##     first, second
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:utils':
## 
##     findMatches
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:data.table':
## 
##     shift
## The following object is masked from 'package:metacoder':
## 
##     reverse
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## The following object is masked from 'package:purrr':
## 
##     reduce
## The following object is masked from 'package:phyloseq':
## 
##     distance
## Loading required package: XVector
## 
## Attaching package: 'XVector'
## The following object is masked from 'package:purrr':
## 
##     compact
## Loading required package: GenomeInfoDb
## 
## Attaching package: 'Biostrings'
## The following object is masked from 'package:metacoder':
## 
##     complement
## The following object is masked from 'package:base':
## 
##     strsplit
## [1] '2.74.1'
library(vegan);packageVersion("vegan")
## Loading required package: permute
## Loading required package: lattice
## [1] '2.6.10'
library(agricolae);packageVersion("agricolae")
## [1] '1.3.7'
library(magrittr);packageVersion("magrittr")
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## [1] '2.0.3'
library(data.table);packageVersion("data.table")
## [1] '1.16.4'
library(DESeq2);packageVersion("DESeq2")
## Loading required package: GenomicRanges
## 
## Attaching package: 'GenomicRanges'
## The following object is masked from 'package:magrittr':
## 
##     subtract
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## The following object is masked from 'package:phyloseq':
## 
##     sampleNames
## [1] '1.46.0'
library(cowplot);packageVersion("cowplot")
## [1] '1.1.3'
library(doParallel);packageVersion("doParallel")
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: iterators
## Loading required package: parallel
## [1] '1.0.17'

Load data from all three experiments

load("Data/ITS/filtered_data/dadsexp.RData")
dads_its.ps<-dads.ps

load("Data/ITS/filtered_data/dads_5asv.RData")
dads_5asv_its.ps<-dads_5asv.ps

load("Data/16S/filtered_data/dadsexp.RData")
dads_16S.ps<-dads.ps

load("Data/16S/filtered_data/dads_5asv.RData")
dads_5asv_16S.ps<-dads_5asv.ps

load("Data/RPS10/filtered_data/dadsexp.RData")
dads_rps10.ps<-dads.ps

load("Data/RPS10/filtered_data/dads_5asv.RData")
dads_5asv_rps10.ps<-dads_5asv.ps

Let’s calculate alpha diversity metrics for ITS dataset

set.seed(1)
variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")

AD_its <- phyloseq::estimate_richness(dads_its.ps, measures = c("Observed", "InvSimpson"))

AD_its <- cbind(AD_its, sample_data(dads_its.ps)[rownames(AD_its), c("Treatment", "Soil", "Incubation")])

variables <- c("Treatment", "Soil", "Incubation")

plots_its <- list()

for (variable in variables) {
  p <- phyloseq::plot_richness(dads_its.ps, x = variable, measures = c("Observed", "InvSimpson")) + 
    geom_boxplot() +
    theme_classic() +  # Use theme_classic() as a starting point
    theme(
      panel.border = element_rect(color = "black", fill = NA, size = 1),  # Add border around plot panel
      axis.line = element_blank(),  # Remove axis lines (they'll be replaced by the panel border)
      plot.margin = unit(c(1, 1, 1, 1), "lines"),  # Add some margin around the plot
      plot.title = element_text(margin = margin(b = 10)),  # Add some space below the title
      axis.title.x = element_text(margin = margin(t = 10)),  # Add some space above the x-axis title
      axis.title.y = element_text(margin = margin(r = 10))  # Add some space to the right of the y-axis title
    )
  plots_its[[variable]] <- p
}
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
  formula <- as.formula(paste(measure, "~", variable))
  aov_result <- aov(formula, data = data)
  tukey_result <- HSD.test(aov_result, variable, group = TRUE)
  letters <- tukey_result$groups$groups
  return(letters)
}

# Perform Tukey HSD tests and store results
tukey_results_ITS <- list()
for (variable in variables) {
  for (measure in measures) {
    tryCatch({
      letters <- get_tukey_letters(AD_its, measure, variable)
      tukey_results_ITS[[paste(measure, variable, sep = "_")]] <- letters
    }, error = function(e) {
      message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
    })
  }
}


print(tukey_results_ITS)
## $Observed_Treatment
## [1] "a" "b"
## 
## $InvSimpson_Treatment
## [1] "a" "b"
## 
## $Observed_Soil
## [1] "a" "a"
## 
## $InvSimpson_Soil
## [1] "a" "a"
## 
## $Observed_Incubation
## [1] "a" "a"
## 
## $InvSimpson_Incubation
## [1] "a" "a"

Let’s calculate alpha diversity metrics for 16S dataset

set.seed(1)
variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")

AD_16S <- phyloseq::estimate_richness(dads_16S.ps, measures = c("Observed", "InvSimpson"))
rownames(AD_16S) <- gsub("\\X", "", rownames(AD_16S))


# Combine richness estimates with sample metadata
AD_16S <- cbind(AD_16S, sample_data(dads_16S.ps)[rownames(AD_16S), c("Treatment", "Soil", "Incubation")])

variables <- c("Treatment", "Soil", "Incubation")

# Create an empty list to store plots
plots_16S <- list()

# Loop through each variable and create the plot
for (variable in variables) {
  p <- phyloseq::plot_richness(dads_16S.ps, x = variable, measures = c("Observed", "InvSimpson")) + 
    geom_boxplot() +
    theme_classic() +  # Use theme_classic() as a starting point
    theme(
      panel.border = element_rect(color = "black", fill = NA, size = 1),  # Add border around plot panel
      axis.line = element_blank(),  # Remove axis lines (they'll be replaced by the panel border)
      plot.margin = unit(c(1, 1, 1, 1), "lines"),  # Add some margin around the plot
      plot.title = element_text(margin = margin(b = 10)),  # Add some space below the title
      axis.title.x = element_text(margin = margin(t = 10)),  # Add some space above the x-axis title
      axis.title.y = element_text(margin = margin(r = 10))  # Add some space to the right of the y-axis title
    )
  plots_16S[[variable]] <- p
}

#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
  formula <- as.formula(paste(measure, "~", variable))
  aov_result <- aov(formula, data = data)
  tukey_result <- HSD.test(aov_result, variable, group = TRUE)
  letters <- tukey_result$groups$groups
  return(letters)
}

# Perform Tukey HSD tests and store results
tukey_results_16S <- list()
for (variable in variables) {
  for (measure in measures) {
    tryCatch({
      letters <- get_tukey_letters(AD_16S, measure, variable)
      tukey_results_16S[[paste(measure, variable, sep = "_")]] <- letters
    }, error = function(e) {
      message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
    })
  }
}


print(tukey_results_16S)
## $Observed_Treatment
## [1] "a" "b"
## 
## $InvSimpson_Treatment
## [1] "a" "b"
## 
## $Observed_Soil
## [1] "a" "a"
## 
## $InvSimpson_Soil
## [1] "a" "a"
## 
## $Observed_Incubation
## [1] "a" "a"
## 
## $InvSimpson_Incubation
## [1] "a" "a"

Let’s calculate alpha diversity metrics for RPS10 dataset

set.seed(1)

variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")

# Estimate alpha diversity
AD_RPS10 <- estimate_richness(dads_rps10.ps, measures = c("Observed", "InvSimpson"))
## Warning in estimate_richness(dads_rps10.ps, measures = c("Observed", "InvSimpson")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
# Combine richness estimates with sample metadata
AD_RPS10 <- cbind(AD_RPS10, sample_data(dads_rps10.ps)[rownames(AD_RPS10), c("Treatment", "Soil", "Incubation")])

variables <- c("Treatment", "Soil", "Incubation")

# Create an empty list to store plots
plots_rps10 <- list()

# Loop through each variable and create the plot
for (variable in variables) {
  p <- phyloseq::plot_richness(dads_rps10.ps, x = variable, measures = c("Observed", "InvSimpson")) + 
    geom_boxplot() +
    theme_classic() +  # Use theme_classic() as a starting point
    theme(
      panel.border = element_rect(color = "black", fill = NA, size = 1),  # Add border around plot panel
      axis.line = element_blank(),  # Remove axis lines (they'll be replaced by the panel border)
      plot.margin = unit(c(1, 1, 1, 1), "lines"),  # Add some margin around the plot
      plot.title = element_text(margin = margin(b = 10)),  # Add some space below the title
      axis.title.x = element_text(margin = margin(t = 10)),  # Add some space above the x-axis title
      axis.title.y = element_text(margin = margin(r = 10))  # Add some space to the right of the y-axis title
    )
  plots_rps10[[variable]] <- p
}
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
  formula <- as.formula(paste(measure, "~", variable))
  aov_result <- aov(formula, data = data)
  tukey_result <- HSD.test(aov_result, variable, group = TRUE)
  letters <- tukey_result$groups$groups
  return(letters)
}

# Perform Tukey HSD tests and store results
tukey_results_RPS10 <- list()
for (variable in variables) {
  for (measure in measures) {
    tryCatch({
      letters <- get_tukey_letters(AD_RPS10, measure, variable)
      tukey_results_RPS10[[paste(measure, variable, sep = "_")]] <- letters
    }, error = function(e) {
      message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
    })
  }
}


print(tukey_results_RPS10)
## $Observed_Treatment
## [1] "a" "a"
## 
## $InvSimpson_Treatment
## [1] "a" "a"
## 
## $Observed_Soil
## [1] "a" "a"
## 
## $InvSimpson_Soil
## [1] "a" "a"
## 
## $Observed_Incubation
## [1] "a" "a"
## 
## $InvSimpson_Incubation
## [1] "a" "a"

Let’s combine the treatment plots

create_label <- function(text, x = 0.5, y = 0.95) {  # Centered the label horizontally
  ggdraw() + 
    draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1)  # Changed hjust to 0.5 for center alignment
}

label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")

# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
  plots_its[["Treatment"]]
)

plots_arrange_16S <- gridExtra::grid.arrange(
  plots_16S[["Treatment"]]
)

plots_arrange_rps10 <- gridExtra::grid.arrange(
  plots_rps10[["Treatment"]]
)

combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width

facetted_figure_trt <- plot_grid(
  combined_plots_16S,
  combined_plots_its,
  combined_plots_rps10,
  nrow = 3,
  align = 'v',
  axis = 'l'
)


facetted_figure_trt_with_margins <- plot_grid(
  NULL, facetted_figure_trt, NULL,
  ncol = 3,
  rel_widths = c(0, 1, 0.08)  # 5% margin on each side
)

facetted_figure_trt_with_margins

# Save the figure
ggsave("Figures/facetted_figure_alphadiv_treatment.svg", facetted_figure_trt_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_treatment.pdf", facetted_figure_trt_with_margins, width=8, height=12)

Let’s combine the soil type plots

create_label <- function(text, x = 0.5, y = 0.95) {  # Centered the label horizontally
  ggdraw() + 
    draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1)  # Changed hjust to 0.5 for center alignment
}

label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")

# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
  plots_its[["Soil"]]
)

plots_arrange_16S <- gridExtra::grid.arrange(
  plots_16S[["Soil"]]
)

plots_arrange_rps10 <- gridExtra::grid.arrange(
  plots_rps10[["Soil"]]
)

combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width

facetted_figure_soil <- plot_grid(
  combined_plots_16S,
  combined_plots_its,
  combined_plots_rps10,
  nrow = 3,
  align = 'v',
  axis = 'l'
)


facetted_figure_soil_with_margins <- plot_grid(
  NULL, facetted_figure_soil, NULL,
  ncol = 3,
  rel_widths = c(0, 1, 0.08)  # 5% margin on each side
)

facetted_figure_soil_with_margins 

# Save the figure
ggsave("Figures/facetted_figure_alphadiv_soil.svg", facetted_figure_soil_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_soil.pdf", facetted_figure_soil_with_margins, width=8, height=12)

Let’s combine the incubation plots

create_label <- function(text, x = 0.5, y = 0.95) {  # Centered the label horizontally
  ggdraw() + 
    draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1)  # Changed hjust to 0.5 for center alignment
}

label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")

# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
  plots_its[["Incubation"]]
)

plots_arrange_16S <- gridExtra::grid.arrange(
  plots_16S[["Incubation"]]
)

plots_arrange_rps10 <- gridExtra::grid.arrange(
  plots_rps10[["Incubation"]]
)

combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1))  # Reduced label width

facetted_figure_inc <- plot_grid(
  combined_plots_16S,
  combined_plots_its,
  combined_plots_rps10,
  nrow = 3,
  align = 'v',
  axis = 'l'
)


facetted_figure_inc_with_margins <- plot_grid(
  NULL, facetted_figure_inc, NULL,
  ncol = 3,
  rel_widths = c(0, 1, 0.08)  # 5% margin on each side
)

facetted_figure_inc_with_margins 

# Save the figure
ggsave("Figures/facetted_figure_alphadiv_incubation.svg", facetted_figure_inc_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_incubation.pdf", facetted_figure_inc_with_margins, width=8, height=12)

Let’s proceed with ordination for the ITS dataset-PCoA

set.seed(1)
strat_rel_dads_its.ps<-transform_sample_counts(dads_5asv_its.ps, function(x) x / sum(x))

ps_PCoA_its = ordinate(strat_rel_dads_its.ps, 
                   method="PCoA", 
                   distance="bray")

plot_its_trt<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_PCoA_its ,
                         color= "Treatment") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_trt<-plot_its_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_its_trt
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

plot_its_inc<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_PCoA_its ,
                         color= "Incubation") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_inc<-plot_its_inc +
   stat_ellipse(linetype =2) +
  theme_bw()

plot_ellipses_its_inc

plot_its_soil<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_PCoA_its ,
                         color= "Soil") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Soil type") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_soil<-plot_its_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_its_soil

Let’s proceed with ordination for the 16S dataset-PCoA

set.seed(1)
strat_rel_dads_16S.ps<-transform_sample_counts(dads_5asv_16S.ps, function(x) x / sum(x))


ps_PCoA_16S = ordinate(strat_rel_dads_16S.ps, 
                   method="PCoA", 
                   distance="bray")

plot_16S_trt<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_PCoA_16S ,
                         color= "Treatment") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_trt<-plot_16S_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_trt

plot_16S_inc<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_PCoA_16S ,
                         color= "Incubation") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_inc<-plot_16S_inc +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_inc

plot_16S_soil<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_PCoA_16S ,
                         color= "Soil") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Soil type") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_soil<-plot_16S_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_soil

Let’s proceed with ordination for the RPS10 dataset-pCoA

set.seed(1)
strat_rel_dads_rps10.ps<-transform_sample_counts(dads_5asv_rps10.ps, function(x) x / sum(x))

ps_PCoA_rps10 = ordinate(strat_rel_dads_rps10.ps, 
                   method="PCoA", 
                   distance="bray")

plot_rps10_trt<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_PCoA_rps10 ,
                         color= "Treatment") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_trt<-plot_rps10_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_trt

plot_rps10_inc<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_PCoA_rps10 ,
                         color= "Incubation") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_inc<-plot_rps10_inc +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_inc

plot_rps10_soil<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_PCoA_rps10 ,
                         color= "Soil") +
  #ggtitle("PCoA ordination plot - relative abundance") +
  labs(color = "Soil type") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_soil<-plot_rps10_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_soil

Let’s combine plots into facetted plot

### Let's combine plots into facetted plot
create_column <- function(ord_plot) {
  plot_grid(
    ord_plot,
    nrow = 1,
    align = 'v',
    axis = 'l',
    rel_heights = c(1)
  )
}

# Create columns for each dataset and factor
col_16S_treatment <- create_column(plot_ellipses_16S_trt)
col_16S_soil <- create_column(plot_ellipses_16S_soil)
col_16S_incubation <- create_column(plot_ellipses_16S_inc)
col_its_treatment <- create_column(plot_ellipses_its_trt)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_its_soil <- create_column(plot_ellipses_its_soil)
col_its_incubation <- create_column(plot_ellipses_its_inc)
col_rps10_treatment <- create_column(plot_ellipses_rps10_trt)
col_rps10_soil <- create_column(plot_ellipses_rps10_soil)
col_rps10_incubation <- create_column(plot_ellipses_rps10_inc)

# Function to create headers with horizontal lines above
create_header <- function(text) {
  ggplot() +
    geom_hline(yintercept = 0.8, color = "black", size = 0.5) +
    geom_text(aes(x = 0.02, y = 0.4, label = text), size = 5, fontface = "bold", hjust = 0) +
    theme_void() +
    theme(plot.margin = margin(5, 0, 0, 10, "pt")) +
    xlim(0, 1) + ylim(0, 1)
}

# Create row plots with labels
row_16S <- plot_grid(col_16S_treatment, col_16S_soil, col_16S_incubation, 
                     ncol = 3, labels = "A.", label_size = 14, hjust = -0.5, vjust = 2)
row_its <- plot_grid(col_its_treatment, col_its_soil, col_its_incubation, 
                     ncol = 3, labels = "B.", label_size = 14, hjust = -0.5, vjust = 2)
row_rps10 <- plot_grid(col_rps10_treatment, col_rps10_soil, col_rps10_incubation, 
                       ncol = 3, labels = "C.", label_size = 14, hjust = -0.5, vjust = 2)

# Combine all elements
facetted_figure <- plot_grid(
  row_16S,
  row_its,
  row_rps10,
  nrow = 3,
  align = 'v',
  axis = 'l',
  rel_heights = c(1.1, 1, 1)  # Slightly more height for the first row
)

facetted_figure

# Save the figure
ggsave("Figures/facetted_figure_with_ordination_PCoA.svg", facetted_figure, width = 12, height = 12)
ggsave("Figures/facetted_figure_with_ordination_PCoA.pdf", facetted_figure, width = 12, height = 12)

Let’s proceed with ordination for the ITS dataset-NMDS

set.seed(1)
strat_rel_dads_its.ps<-transform_sample_counts(dads_5asv_its.ps, function(x) x / sum(x))

ps_NMDS_its = ordinate(strat_rel_dads_its.ps, 
                   method="NMDS", 
                   distance="bray")
## Run 0 stress 9.209135e-05 
## Run 1 stress 9.711561e-05 
## ... Procrustes: rmse 4.557734e-05  max resid 0.000145828 
## ... Similar to previous best
## Run 2 stress 9.941333e-05 
## ... Procrustes: rmse 4.783054e-05  max resid 0.0001492673 
## ... Similar to previous best
## Run 3 stress 9.765011e-05 
## ... Procrustes: rmse 5.327067e-05  max resid 0.0001430395 
## ... Similar to previous best
## Run 4 stress 9.298127e-05 
## ... Procrustes: rmse 9.231115e-05  max resid 0.0002031309 
## ... Similar to previous best
## Run 5 stress 9.956743e-05 
## ... Procrustes: rmse 6.453389e-05  max resid 0.0001518031 
## ... Similar to previous best
## Run 6 stress 9.81565e-05 
## ... Procrustes: rmse 6.29341e-05  max resid 0.0001313132 
## ... Similar to previous best
## Run 7 stress 9.415583e-05 
## ... Procrustes: rmse 2.284447e-05  max resid 8.763404e-05 
## ... Similar to previous best
## Run 8 stress 9.704892e-05 
## ... Procrustes: rmse 0.0001012395  max resid 0.0002088867 
## ... Similar to previous best
## Run 9 stress 9.824374e-05 
## ... Procrustes: rmse 6.666882e-05  max resid 0.0001454821 
## ... Similar to previous best
## Run 10 stress 9.985598e-05 
## ... Procrustes: rmse 0.000100697  max resid 0.0002126261 
## ... Similar to previous best
## Run 11 stress 9.996398e-05 
## ... Procrustes: rmse 6.547809e-05  max resid 0.0001601838 
## ... Similar to previous best
## Run 12 stress 9.761524e-05 
## ... Procrustes: rmse 4.580355e-05  max resid 0.0001486625 
## ... Similar to previous best
## Run 13 stress 9.698329e-05 
## ... Procrustes: rmse 5.20759e-05  max resid 0.0001371574 
## ... Similar to previous best
## Run 14 stress 9.946996e-05 
## ... Procrustes: rmse 5.137318e-05  max resid 0.0001439045 
## ... Similar to previous best
## Run 15 stress 0.0002273924 
## ... Procrustes: rmse 0.000322125  max resid 0.001292199 
## ... Similar to previous best
## Run 16 stress 9.993753e-05 
## ... Procrustes: rmse 4.747406e-05  max resid 0.0001496008 
## ... Similar to previous best
## Run 17 stress 9.9414e-05 
## ... Procrustes: rmse 7.479917e-05  max resid 0.0001668549 
## ... Similar to previous best
## Run 18 stress 9.94672e-05 
## ... Procrustes: rmse 0.00010456  max resid 0.0002113862 
## ... Similar to previous best
## Run 19 stress 9.748534e-05 
## ... Procrustes: rmse 6.363321e-05  max resid 0.0001336849 
## ... Similar to previous best
## Run 20 stress 9.760704e-05 
## ... Procrustes: rmse 6.106859e-05  max resid 0.0001475123 
## ... Similar to previous best
## *** Best solution repeated 20 times
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_its_trt<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_NMDS_its ,
                         color= "Treatment") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_trt<-plot_its_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_its_trt

plot_its_inc<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_NMDS_its ,
                         color= "Incubation") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_inc<-plot_its_inc +
   stat_ellipse(linetype =2) +
  theme_bw()

plot_ellipses_its_inc

plot_its_soil<- plot_ordination(strat_rel_dads_its.ps, 
                         ps_NMDS_its ,
                         color= "Soil") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Soil type") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_its_soil<-plot_its_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_its_soil

Let’s proceed with ordination for the 16S dataset-NMDS

set.seed(1)
strat_rel_dads_16S.ps<-transform_sample_counts(dads_5asv_16S.ps, function(x) x / sum(x))


ps_NMDS_16S = ordinate(strat_rel_dads_16S.ps, 
                   method="NMDS", 
                   distance="bray")
## Run 0 stress 8.996366e-05 
## Run 1 stress 8.317019e-05 
## ... New best solution
## ... Procrustes: rmse 0.003178651  max resid 0.01480472 
## Run 2 stress 9.962288e-05 
## ... Procrustes: rmse 0.02913161  max resid 0.1366179 
## Run 3 stress 9.148e-05 
## ... Procrustes: rmse 0.03153671  max resid 0.1479125 
## Run 4 stress 9.094561e-05 
## ... Procrustes: rmse 0.02732531  max resid 0.1281317 
## Run 5 stress 8.507797e-05 
## ... Procrustes: rmse 0.03347368  max resid 0.1570019 
## Run 6 stress 9.458946e-05 
## ... Procrustes: rmse 0.01482677  max resid 0.06938939 
## Run 7 stress 9.602634e-05 
## ... Procrustes: rmse 0.03118334  max resid 0.1427293 
## Run 8 stress 9.124274e-05 
## ... Procrustes: rmse 0.008086828  max resid 0.03753532 
## Run 9 stress 9.099226e-05 
## ... Procrustes: rmse 0.007628467  max resid 0.03563048 
## Run 10 stress 9.193224e-05 
## ... Procrustes: rmse 0.03255636  max resid 0.152698 
## Run 11 stress 8.232243e-05 
## ... New best solution
## ... Procrustes: rmse 0.04187255  max resid 0.1900619 
## Run 12 stress 9.230235e-05 
## ... Procrustes: rmse 0.06475307  max resid 0.303561 
## Run 13 stress 8.965179e-05 
## ... Procrustes: rmse 0.0480184  max resid 0.2242121 
## Run 14 stress 9.391617e-05 
## ... Procrustes: rmse 0.06917108  max resid 0.3244004 
## Run 15 stress 9.920624e-05 
## ... Procrustes: rmse 0.0691793  max resid 0.3244382 
## Run 16 stress 9.5959e-05 
## ... Procrustes: rmse 0.06341095  max resid 0.2972143 
## Run 17 stress 8.293653e-05 
## ... Procrustes: rmse 0.06930239  max resid 0.3250205 
## Run 18 stress 8.42273e-05 
## ... Procrustes: rmse 0.06202637  max resid 0.2906629 
## Run 19 stress 8.432815e-05 
## ... Procrustes: rmse 0.06944932  max resid 0.32571 
## Run 20 stress 9.080799e-05 
## ... Procrustes: rmse 0.02793296  max resid 0.1292364 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_16S_trt<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_NMDS_16S ,
                         color= "Treatment") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_trt<-plot_16S_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_trt
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

plot_16S_inc<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_NMDS_16S ,
                         color= "Incubation") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_inc<-plot_16S_inc +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_inc
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

plot_16S_soil<- plot_ordination(strat_rel_dads_16S.ps, 
                         ps_NMDS_16S ,
                         color= "Soil") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Soil type") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_16S_soil<-plot_16S_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_16S_soil
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

Let’s proceed with ordination for the RPS10 dataset-NMDS

set.seed(1)
strat_rel_dads_rps10.ps<-transform_sample_counts(dads_5asv_rps10.ps, function(x) x / sum(x))

ps_NMDS_rps10 = ordinate(strat_rel_dads_rps10.ps, 
                   method="NMDS", 
                   distance="bray")
## Run 0 stress 9.646111e-05 
## Run 1 stress 9.681011e-05 
## ... Procrustes: rmse 0.0001489952  max resid 0.0004893922 
## ... Similar to previous best
## Run 2 stress 8.725043e-05 
## ... New best solution
## ... Procrustes: rmse 9.343795e-05  max resid 0.0003534458 
## ... Similar to previous best
## Run 3 stress 9.218624e-05 
## ... Procrustes: rmse 9.166428e-05  max resid 0.0002957858 
## ... Similar to previous best
## Run 4 stress 9.131775e-05 
## ... Procrustes: rmse 8.62371e-05  max resid 0.0003601615 
## ... Similar to previous best
## Run 5 stress 9.687112e-05 
## ... Procrustes: rmse 8.700801e-05  max resid 0.0001959933 
## ... Similar to previous best
## Run 6 stress 9.390856e-05 
## ... Procrustes: rmse 0.0001078327  max resid 0.0003442075 
## ... Similar to previous best
## Run 7 stress 9.920391e-05 
## ... Procrustes: rmse 5.200916e-05  max resid 0.000175748 
## ... Similar to previous best
## Run 8 stress 9.720681e-05 
## ... Procrustes: rmse 9.55066e-05  max resid 0.0002879497 
## ... Similar to previous best
## Run 9 stress 8.898622e-05 
## ... Procrustes: rmse 9.951896e-05  max resid 0.0002879941 
## ... Similar to previous best
## Run 10 stress 8.54193e-05 
## ... New best solution
## ... Procrustes: rmse 7.24128e-05  max resid 0.0002186124 
## ... Similar to previous best
## Run 11 stress 9.003123e-05 
## ... Procrustes: rmse 6.651538e-05  max resid 0.0002296343 
## ... Similar to previous best
## Run 12 stress 8.8711e-05 
## ... Procrustes: rmse 6.639742e-05  max resid 0.0002112206 
## ... Similar to previous best
## Run 13 stress 9.142818e-05 
## ... Procrustes: rmse 0.0001033687  max resid 0.0002963042 
## ... Similar to previous best
## Run 14 stress 8.759001e-05 
## ... Procrustes: rmse 0.0001045275  max resid 0.0002975806 
## ... Similar to previous best
## Run 15 stress 8.993688e-05 
## ... Procrustes: rmse 9.302378e-05  max resid 0.0002955304 
## ... Similar to previous best
## Run 16 stress 8.753563e-05 
## ... Procrustes: rmse 7.733486e-05  max resid 0.0002073648 
## ... Similar to previous best
## Run 17 stress 8.508021e-05 
## ... New best solution
## ... Procrustes: rmse 7.811403e-05  max resid 0.0002102992 
## ... Similar to previous best
## Run 18 stress 9.858152e-05 
## ... Procrustes: rmse 9.462793e-05  max resid 0.0002329584 
## ... Similar to previous best
## Run 19 stress 9.620554e-05 
## ... Procrustes: rmse 9.740211e-05  max resid 0.000223294 
## ... Similar to previous best
## Run 20 stress 9.002458e-05 
## ... Procrustes: rmse 6.780994e-05  max resid 0.0001783893 
## ... Similar to previous best
## *** Best solution repeated 4 times
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_rps10_trt<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_NMDS_rps10 ,
                         color= "Treatment") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Treatment") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_trt<-plot_rps10_trt +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_trt

plot_rps10_inc<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_NMDS_rps10 ,
                         color= "Incubation") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Incubation") + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_inc<-plot_rps10_inc +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_inc

plot_rps10_soil<- plot_ordination(strat_rel_dads_rps10.ps, 
                         ps_NMDS_rps10 ,
                         color= "Soil") +
  #ggtitle("NMDS ordination plot - relative abundance") +
  labs(color = "Soil type") +
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title =  element_text(size=13), legend.text = element_text(size = 18))

plot_ellipses_rps10_soil<-plot_rps10_soil +
   stat_ellipse(linetype = 2) +
  theme_bw()

plot_ellipses_rps10_soil

Let’s combine plots into facetted plot

create_column <- function(ord_plot) {
  plot_grid(
    ord_plot,
    nrow = 1,
    align = 'v',
    axis = 'l',
    rel_heights = c(1)
  )
}

# Create columns for each dataset and factor
col_16S_treatment <- create_column(plot_ellipses_16S_trt)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_16S_soil <- create_column(plot_ellipses_16S_soil)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_16S_incubation <- create_column(plot_ellipses_16S_inc)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_its_treatment <- create_column(plot_ellipses_its_trt)
col_its_soil <- create_column(plot_ellipses_its_soil)
col_its_incubation <- create_column(plot_ellipses_its_inc)
col_rps10_treatment <- create_column(plot_ellipses_rps10_trt)
col_rps10_soil <- create_column(plot_ellipses_rps10_soil)
col_rps10_incubation <- create_column(plot_ellipses_rps10_inc)

# Create row plots with labels
row_16S <- plot_grid(col_16S_treatment, col_16S_soil, col_16S_incubation, 
                     ncol = 3, labels = "A.", label_size = 14, hjust = -0.5, vjust = 2)
row_its <- plot_grid(col_its_treatment, col_its_soil, col_its_incubation, 
                     ncol = 3, labels = "B.", label_size = 14, hjust = -0.5, vjust = 2)
row_rps10 <- plot_grid(col_rps10_treatment, col_rps10_soil, col_rps10_incubation, 
                       ncol = 3, labels = "C.", label_size = 14, hjust = -0.5, vjust = 2)

# Combine all elements
facetted_figure <- plot_grid(
  row_16S,
  row_its,
  row_rps10,
  nrow = 3,
  align = 'v',
  axis = 'l',
  rel_heights = c(1.1, 1, 1)  # Slightly more height for the first row
)

facetted_figure

# Save the figure
ggsave("Figures/facetted_figure_with_ordination_NMDS.svg", facetted_figure, width = 12, height = 12)
ggsave("Figures/facetted_figure_with_ordination_NMDS.pdf", facetted_figure, width = 12, height = 12)

Software used

sessioninfo::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.2 (2024-10-31)
##  os       macOS Ventura 13.2.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/Los_Angeles
##  date     2025-02-19
##  pandoc   3.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
##  quarto   1.3.450 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version  date (UTC) lib source
##  abind                  1.4-8    2024-09-12 [1] CRAN (R 4.4.1)
##  ade4                   1.7-23   2025-02-14 [1] CRAN (R 4.4.1)
##  agricolae            * 1.3-7    2023-10-22 [1] CRAN (R 4.4.0)
##  AlgDesign              1.2.1.1  2024-09-21 [1] CRAN (R 4.4.1)
##  ape                    5.8-1    2024-12-16 [1] CRAN (R 4.4.1)
##  Biobase              * 2.66.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocGenerics         * 0.52.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocParallel           1.40.0   2024-11-23 [1] Bioconductor 3.20 (R 4.4.2)
##  biomformat             1.34.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  Biostrings           * 2.74.1   2024-12-16 [1] Bioconductor 3.20 (R 4.4.2)
##  bslib                  0.9.0    2025-01-30 [1] CRAN (R 4.4.1)
##  cachem                 1.1.0    2024-05-16 [1] CRAN (R 4.4.1)
##  cli                    3.6.4    2025-02-13 [1] CRAN (R 4.4.1)
##  cluster                2.1.8    2024-12-11 [1] CRAN (R 4.4.1)
##  codetools              0.2-20   2024-03-31 [1] CRAN (R 4.4.2)
##  colorspace             2.1-1    2024-07-26 [1] CRAN (R 4.4.1)
##  cowplot              * 1.1.3    2024-01-22 [1] CRAN (R 4.4.0)
##  crayon                 1.5.3    2024-06-20 [1] CRAN (R 4.4.1)
##  data.table           * 1.16.4   2024-12-06 [1] CRAN (R 4.4.1)
##  DelayedArray           0.32.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  DESeq2               * 1.46.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  digest                 0.6.37   2024-08-19 [1] CRAN (R 4.4.1)
##  doParallel           * 1.0.17   2022-02-07 [1] CRAN (R 4.4.0)
##  dplyr                * 1.1.4    2023-11-17 [1] CRAN (R 4.4.0)
##  evaluate               1.0.3    2025-01-10 [1] CRAN (R 4.4.1)
##  farver                 2.1.2    2024-05-13 [1] CRAN (R 4.4.1)
##  fastmap                1.2.0    2024-05-15 [1] CRAN (R 4.4.1)
##  foreach              * 1.5.2    2022-02-02 [1] CRAN (R 4.4.0)
##  furrr                * 0.3.1    2022-08-15 [1] CRAN (R 4.4.0)
##  future               * 1.34.0   2024-07-29 [1] CRAN (R 4.4.0)
##  generics               0.1.3    2022-07-05 [1] CRAN (R 4.4.1)
##  GenomeInfoDb         * 1.42.3   2025-01-27 [1] Bioconductor 3.20 (R 4.4.2)
##  GenomeInfoDbData       1.2.13   2025-02-07 [1] Bioconductor
##  GenomicRanges        * 1.58.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  ggplot2              * 3.5.1    2024-04-23 [1] CRAN (R 4.4.0)
##  globals                0.16.3   2024-03-08 [1] CRAN (R 4.4.1)
##  glue                   1.8.0    2024-09-30 [1] CRAN (R 4.4.1)
##  gridExtra              2.3      2017-09-09 [1] CRAN (R 4.4.1)
##  gtable                 0.3.6    2024-10-25 [1] CRAN (R 4.4.1)
##  hms                    1.1.3    2023-03-21 [1] CRAN (R 4.4.0)
##  htmltools              0.5.8.1  2024-04-04 [1] CRAN (R 4.4.1)
##  httr                   1.4.7    2023-08-15 [1] CRAN (R 4.4.0)
##  igraph                 2.1.4    2025-01-23 [1] CRAN (R 4.4.1)
##  IRanges              * 2.40.1   2024-12-05 [1] Bioconductor 3.20 (R 4.4.2)
##  iterators            * 1.0.14   2022-02-05 [1] CRAN (R 4.4.1)
##  jquerylib              0.1.4    2021-04-26 [1] CRAN (R 4.4.0)
##  jsonlite               1.8.9    2024-09-20 [1] CRAN (R 4.4.1)
##  knitr                  1.49     2024-11-08 [1] CRAN (R 4.4.1)
##  labeling               0.4.3    2023-08-29 [1] CRAN (R 4.4.1)
##  lattice              * 0.22-6   2024-03-20 [1] CRAN (R 4.4.2)
##  lifecycle              1.0.4    2023-11-07 [1] CRAN (R 4.4.1)
##  listenv                0.9.1    2024-01-29 [1] CRAN (R 4.4.1)
##  locfit                 1.5-9.11 2025-02-03 [1] CRAN (R 4.4.1)
##  magrittr             * 2.0.3    2022-03-30 [1] CRAN (R 4.4.1)
##  MASS                   7.3-64   2025-01-04 [1] CRAN (R 4.4.1)
##  Matrix                 1.7-2    2025-01-23 [1] CRAN (R 4.4.1)
##  MatrixGenerics       * 1.18.1   2025-01-09 [1] Bioconductor 3.20 (R 4.4.2)
##  matrixStats          * 1.5.0    2025-01-07 [1] CRAN (R 4.4.1)
##  metacoder            * 0.3.8    2025-02-11 [1] CRAN (R 4.4.1)
##  mgcv                   1.9-1    2023-12-21 [1] CRAN (R 4.4.2)
##  multtest               2.62.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  munsell                0.5.1    2024-04-01 [1] CRAN (R 4.4.1)
##  nlme                   3.1-167  2025-01-27 [1] CRAN (R 4.4.1)
##  parallelly             1.42.0   2025-01-30 [1] CRAN (R 4.4.1)
##  permute              * 0.9-7    2022-01-27 [1] CRAN (R 4.4.1)
##  phyloseq             * 1.50.0   2024-10-29 [1] Bioconductor 3.20 (R 4.4.2)
##  pillar                 1.10.1   2025-01-07 [1] CRAN (R 4.4.1)
##  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 4.4.1)
##  plyr                   1.8.9    2023-10-02 [1] CRAN (R 4.4.1)
##  purrr                * 1.0.4    2025-02-05 [1] CRAN (R 4.4.1)
##  R6                     2.6.1    2025-02-15 [1] CRAN (R 4.4.1)
##  ragg                   1.3.3    2024-09-11 [1] CRAN (R 4.4.1)
##  Rcpp                   1.0.14   2025-01-12 [1] CRAN (R 4.4.1)
##  readr                * 2.1.5    2024-01-10 [1] CRAN (R 4.4.0)
##  reshape2               1.4.4    2020-04-09 [1] CRAN (R 4.4.0)
##  rhdf5                  2.50.2   2025-01-09 [1] Bioconductor 3.20 (R 4.4.2)
##  rhdf5filters           1.18.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  Rhdf5lib               1.28.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  rlang                  1.1.5    2025-01-17 [1] CRAN (R 4.4.1)
##  rmarkdown              2.29     2024-11-04 [1] CRAN (R 4.4.1)
##  rstudioapi             0.17.1   2024-10-22 [1] CRAN (R 4.4.1)
##  S4Arrays               1.6.0    2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  S4Vectors            * 0.44.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  sass                   0.4.9    2024-03-15 [1] CRAN (R 4.4.0)
##  scales                 1.3.0    2023-11-28 [1] CRAN (R 4.4.0)
##  sessioninfo            1.2.3    2025-02-05 [1] CRAN (R 4.4.1)
##  SparseArray            1.6.1    2025-01-20 [1] Bioconductor 3.20 (R 4.4.2)
##  stringi                1.8.4    2024-05-06 [1] CRAN (R 4.4.1)
##  stringr              * 1.5.1    2023-11-14 [1] CRAN (R 4.4.0)
##  SummarizedExperiment * 1.36.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  survival               3.8-3    2024-12-17 [1] CRAN (R 4.4.1)
##  svglite                2.1.3    2023-12-08 [1] CRAN (R 4.4.0)
##  systemfonts            1.2.1    2025-01-20 [1] CRAN (R 4.4.1)
##  textshaping            1.0.0    2025-01-20 [1] CRAN (R 4.4.1)
##  tibble                 3.2.1    2023-03-20 [1] CRAN (R 4.4.0)
##  tidyselect             1.2.1    2024-03-11 [1] CRAN (R 4.4.0)
##  tzdb                   0.4.0    2023-05-12 [1] CRAN (R 4.4.0)
##  UCSC.utils             1.2.0    2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  vctrs                  0.6.5    2023-12-01 [1] CRAN (R 4.4.0)
##  vegan                * 2.6-10   2025-01-29 [1] CRAN (R 4.4.1)
##  withr                  3.0.2    2024-10-28 [1] CRAN (R 4.4.1)
##  xfun                   0.50     2025-01-07 [1] CRAN (R 4.4.1)
##  XVector              * 0.46.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##  yaml                   2.3.10   2024-07-26 [1] CRAN (R 4.4.1)
##  zlibbioc               1.52.0   2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
##  * ── Packages attached to the search path.
## 
## ──────────────────────────────────────────────────────────────────────────────